'''
This script filter out dummy sequences with low entropy
'''

import os
import getopt
import sys
from parse_fasta_file import get_tag_seq_from_fasta_w_substr

def entropy( seq ):
    '''Currently return how many consecutive rep di-nt is 
       TODO: replace this by returning the true entropy  
    '''
    counter = 0
    for i in range( len( seq ) - 1 ):
        if seq[i] == seq[i + 1]:
            counter += 1
    return counter

def filter_reads( iFileN, oFileN, threshold ):
    iFile = open( iFileN, "r" )
    oFile = open( "removed" + oFileN, "w" )
    oFile2 = open( oFileN, "w" )

    dummyReadCount = 0
    nonDummyReadCount = 0
    while True:
        [seqN, seq] = ["", ""]
        [seqN, seq] = get_tag_seq_from_fasta_w_substr( iFile, '>' )
        # print seqN
        if seqN:
            dummyScore = entropy( seq )
            if dummyScore < threshold:
                oFile.write( '>' + seqN + '\n' )
                oFile.write( seq + '\n' )
                nonDummyReadCount += 1
            else:
                oFile2.write( '>' + seqN + ' ' + str( dummyScore ) + '\n' )
                oFile2.write( seq + '\n' )
                dummyReadCount += 1
        else:
            line = iFile.readline()
            if line:
                print "Format in ", line, " is not expected."
            break
    iFile.close()
    oFile.close()
    oFile2.close()
    print "There are ", dummyReadCount, "dummy reads."
    print "There are ", nonDummyReadCount, "non dummy reads."

def usage():
    sys.stderr.write( "Working dir is " + os.getcwd() + '\n' )
    print "usage: ", sys.argv[0]
    print "The two parameters should be "
    print "(1) the input .fa file path"
    print "(2) the output .fa file path"
    print "(3) the threshold for filtering"
    sys.exit( 2 )

def main():
    try:
        opts, args = getopt.getopt( sys.argv[1:], 'h', ["help"] )
    except getopt.error:
        usage()
        sys.exit( 2 );
    # process options and argument 
    for o in opts:
        if o in ( "-h", "--help" ):
            usage()
            sys.exit();
    if len( args ) == 3:
        iFileN = args[0]
        oFileN = args[1]
        threshold = int( args[2] )
        filter_reads( iFileN, oFileN, threshold )
    else:
        usage()

if __name__ == '__main__':
    main()
